Consistent Velocity Approximation for FiniteVolume or Element Discretizations of DensityDriven Flow in Porous Media

Many controversial numerical simulations of density driven ow problems presented in di erent papers suggest the necessity of careful and proper design of discretization. One of the possible numerical artifact that appears if standard nite element or nite volume methods are formally applied, is the presence of arti cial numerical velocities even for the case of hydrostatic conditions. This subject is known as a \consistent velocity approximation" Voss and Souza [1987] where a restrictive assumption on the interpolation of pressure and density must be used in the computations of velocity from (1) ~v = K (rp (c)~g): An algorithm that is so far used for the determination of consistent velocity, can be validated only for a linearly variable density and only in the case of one{ dimensional hydrostatic equilibrium. We present a new formulation of consistent velocity approximation for the arbitrary choice of discretization technique. Such formulation can be validated exactly for a general nonlinear dependence of density with no assumptions about hydrostatic conditions. Further, a simple relation to previous consistency formulation is given with the conclusion that it is equivalent for the linearly variable density. This enables us to present the algorithm also for nite elements not found in literature (triangles, tetrahedras, : : :). The new idea of consistent velocity approximation was used in the context of nite volume discretization. Together with an upwind scheme it results in stable, physically reasonable and grid nondependent behaviour of numerical solutions. 1. Mathematical model We want to solve for the unknown functions c = c(x; t) and p = p(x; t) considered on (0; T ) with R2 or 3 and (0; T ) R the following problem: (1.1) @t( c) +r ( c~v D rc) = Qc; (1.2) @t +r ( ~v) = Qp where the uid velocity ~v is explicitly given by (1.3) ~v = K (rp ~g): Standard boundary and initial conditions are supposed. The variable density is expressed through a functional dependance on the concentration (1.4) = (c): The additional nonlinear coupling of (1.1) and (1.2) is introduced through the dependence of viscosity = (c) on the concentration and the dependence of dispersion part in the di usion{dispersion tensor D = D(~v) on the velocity. Moreover the di usion{dispersion tensor D, the permeability tensor K and the porosity are allowed to vary in space. By ~g we denote the gravity vector. Although the source terms Qc and Qp should be studied in a very general form, we simplify our treatment by allowing the dependence only on space. As this model is extensively studied in literature we do not present any further details, ignoring the fact of their importance. 2. Consistent velocity approximation Before applying a nite volume discretization to solve numerically the equations (1.1){(1.2), we want to clear out a problem that is known from literature as a \consistent velocity aproximation". The most typical explanation of the topic is through the presence of arti cial numerical velocities for the hydrostatic conditions if not a consistent representation for pressure and density is used. We demonstrate it on a simple example where the density is varying linearly in the direction of gravity (2.1) = (z) = 0 + ( 1 0)z; z 2 [0; 1] where the interval [0; 1] was choosen for a simplicity. To reach the zero velocity in (1.3) we need (2.2) @xp = @yp = 0; @zp = g where g is the gravity constant. From the last equation in (2.2) we get (2.3) p = p(z) = p0 g 0z + 1 0 2 z2 : For z = 1 we denote p1 := p(1) = p0 g ( 0 + 1)=2 . Clearly for the case of the linear density (2.1) to maintain a hydrostatic equilibrium the pressure 1 should exhibit the quadratic shape (2.3). A common assumption for many discretization algorithms is that the approximated pressure ~ p is linearly interpolated between two given values, so for the exact nodal values p0 and p1 we have (2.4) ~ p = p0 + (p1 p0)z; z 2 [0; 1]. The computations of velocities using (2.4) in (1.3) lead to in general nonzero z{component of approximated velocity: (2.5) ~ vz = kzz g( 0 1)(1=2 z); z 2 [0; 1]: For a nonconstant linear density the velocity vanishes only at the middle point of the interval and takes the maximum value with opposite signs at the left and right point. Using such nonconsistent velocities can result in an arti cial salt distribution through the convective term and the velocity{dependent dispersion tensor in (1.2) and consequently can in uence the ow itself through the changes of density in (1.1). To reduce this numerical artifact to a minimum one has to use a very ne discretization length in z direction. The remedy of nonconsistent velocity approximation is well documented in the work of Voss and Souza [1987] including the discussion of possible solutions in the framework of nite element methods. One of the conclusion was that either quadratic shape functions shall be used for a pressure approximation or the representation of density shall be reduced to a piecewise constant function. The rst idea results in additional computational costs and the second one in the precision reduction. For the case of the nite element method on a quadrilateral mesh with bilinear ansatz functions Voss and Souza [1987] suggest a compromis procedure. As we present this algorithm comprehensively in the next section, here we introduce it only very brie y. The authors suggest to compute \a local discretisation of ~g with a spatial functionality that is consistent with local presure derivatives". Using then the local representation by means of the reference element [0; 1] [0; 1]; see the next section, their formulas take the form: (2.6) ( ~g) ( ) = 1 2 4 P i=1 ig i @Ni @ , ( ~g) ( ) = 1 2 4 P i=1 ig i @Ni @ : Here i and ~gi are the values of density and of gravity in the local coordinates ( ; ) at the nodes of an element and Ni = Ni( ; ) are the corresponding local ansatz functions, see the next section. The authors claim that the validity of (2.6) can be found by inspection and that no particular signi cance should be attached to the absolute values of basis function derivatives in (2.6). Leijnse [1992] presents a more suitable and equivalent form of (2.6): 2 (2.7) ( ~g) ( ) = 12( 4 P i=1 ig iNi(0; ) + 4 P i=1 ig iNi(1; )); (2.8) ( ~g) ( ) = 12( 4 P i=1 ig iNi( ; 0) + 4 P i=1 ig iNi( ; 1)): Now the algorithm can be understood as an averaging of density{gravity term in apropriate direction to insure the same \spatial functionality" as for the gradient of pressure. Nevertheless there are several limitations of this approach. First these formulas are given ad hoc and in fact (2.6) is presented in Voss and Souza [1987] without 1/2 as they use the reference element [ 1; 1] [ 1; 1] and the appearance of (2.6) is di erent with the di erent representation of ansatz functions. The derivation from hydrostatic conditions in a two dimensional case fails for a general density = (x; z) as here there is no hydrostatic equilibrium due to the nonexistence of potential in the presence of variable density{gravity term. The validation of (2.6) can be provided only for the case of linearly variable density and a generalization of (2.6) is not clear for a nonlinear form of (1.4) that is often used, see Leijnse [1992]. All these shortcomings we overcome by a new formulation of consistent velocity in the next section. 3. New formulation of consistent velocity We restrict here ourselves to the two{dimensional case as a generalization forR3 is straightforward. We suppose that the domain is triangulated by the nite elements that are isoparametrically equivalent to a reference element with local coordinates ( ; ). We consider a special form of the mapping from the local coordinates ( ; ) to the global ones (x; z) (3.1) x = x( ; ) = K P k=1xkNk( ; ), z = z( ; ) = K P k=1 zkNk( ; ) where (xk; zk), k = 1; : : : ; K are the coordinates of all vertices of the element considered and Nk are the local ansatz functions. Then by denoting J to be a Jacobian matrix (3.2) J = J( ; ) = K P k=1 xk@ Nk zk@ Nk xk@ Nk zk@ Nk ! we have the equivalent formulation of (1.3) in local coordinates (3.3) ~v = ~v(x( ; ); z( ; )) = K J 1 (r( ; )p J ~g). Finally considering the notation (3.4) g g ! := g ( ; ) g ( ; ) ! = J ~g we introduce (3.5) h = h ( ; ) = R0 (s; )g (s; )ds, 3 (3.6) h = h ( ; ) = R0 ( ; s)g ( ; s)ds. The motivation of (3.5){(3.6) is clear from the equivalent form of (3.3) (3.7) ~v = K J 1 @ (p h ) @ (p h ) !. With (3.7) the consistent approximation of ~v arises naturally for any discretization. For given or computed values of concentration at the vertices of an element together with a speci c (polynomial) ansatz, the values of the functions h have to be determined according to (1.4) and (3.5){(3.6) and the same (polynomial) ansatz as for the pressure has to be applied then. We present details for the two most common type of elements and interpolations in the 2D case. We start with the reference triangle having the corners (0; 0) (1; 0) and (0; 1), using the local linear ansatz functions: (3.8) N1( ; ) = 1 ; N2( ; ) = ; N3( ; ) = : The Jacobian J from (3.2) does not depend then on the local coordinates and the gravity in local coordinates from (3.4) is a constant vector. By ck; k = 1; : : : ; K we denote the values of c at the vertices (xk; zk) and we suppose the interpolation of the form (3.9) c = c(x( ; ); z( ; )) = K P k=1 ckNk( ; ): Next we compute the values of h and h using (3.5){(3.6) at all vertices of the element denoting them by h k and h k. We get rst: (3.10) h 1 = h 3 = 0, h 2 = g 1 R0 (c(s; 0))ds. Further from (3.8) and (3.9) we have (3.11) 1 R0 (c(s; 0))ds = 1 R0 (c1 + (c2 c1)s)ds = 1 c2 c1 c2 R c1 (s)ds: If we denote by R the primitive of : (3.12) R = R(c) = R (c)dc then we end with (3.13) h 2 = g R(c2) R(c1) c2 c1 . Similarly (3.14) h 1 = h 2 = 0, h 3 = g 1 R0 (c(0; s))ds = g R(c3) R(c1) c3 c1 . If next for the pressure the interpolation is supposed of the form (3.15) p = p(x( ; ); z( ; )) = K P k=1 pkNk( ; ) then the same ansatz for h and h has to be used. First we introduce the new formulation of consistent velocity approximation in a general form: 4 (3.16) ~v(x( ; ); z( ; )) = K J 1 K P k=1 (pk h k)@ Nk( ; ) (pk h k)@ Nk( ; ) !. To compare it with the \standard" local discretization of density{gravity term (2.6) of Voss and Souza [1987] and (2.7){(2.8) of Leijnse [1992], we can consider the density-gravity term ~g in the local coordinates to be given by (3.17) ( ~g) ( ; ) = @ h ( ; ) = K P k=1h k@ Nk( ; ), (3.18) ( ~g) ( ; ) = @ h ( ; ) = K P k=1h k@ Nk( ; ). For the reference the triangle we have that @ h h 2 and @ h h 3 so the formulas (3.17){(3.18) take the simple form (3.19) ( ~g) ( ; ) h 2 = R(c2) R(c1) c2 c1 g , ( ~g) ( ; ) h 3 = g R(c3) R(c1) c3 c1 . The relation to the \standard" consistent velocity formulation can be now introduced. By considering the numerical integration in (3.10) and (3.14) by the trapezoidal rule, we get: (3.20) h 2 1 2( (c1) + (c2))g , h 3 1 2( (c1) + (c3))g . The formulas (3.20) together with (3.16) can be considered as the consistent velocity approximation in the traditional sense for triangle element that we could not nd in literature. For the case of a linear dependance of the density on concentration it is equivalent, but in general only approximative. This interpretation is fully con rmed with the known results for the case of quadrilateral element. Here the reference element have the corners (0; 0), (1; 0), (1; 1) and (0; 1) with the bilinear local ansatz functions: (3.21) N1 = (1 )(1 ), N2 = (1 ), N3 = , N4 = (1 ) . The Jacobian J in (3.2) is in general space dependant but the gravity in the local coordinates from (3.4) takes a special form: (3.22) g g ! = g ( ) g ( ) ! : By computing the values of h and h at the corners we obtain (3.23) h 1 = h 4 = 0; h 2 = g (0) 1 R0 (c(s; 0))ds = g (0)R(c2) R(c1) c2 c1 h 3 = g (1) 1 R0 (c(s; 1))ds = g (1)R(c3) R(c4) c3 c4 (3.24) h 1 = h 2 = 0; h 3 = g (1) 1 R0 (c(1; s))ds = g (1)R(c3) R(c2) c3 c2 h 4 = g (0) 1 R0 (c(0; s))ds = g (0)R(c4) R(c1) c4 c1 The new local discretization of density{gravity takes the form (3.25) ( ~g) = @ h = h 2(1 ) + h 3 , ( ~g) = @ h = h 3 + h 4(1 ). Again using the numerical integration in (3.23) and (3.24) we get 5 (3.26) h 212( (c1) + (c2))g (0); h 312( (c3) + (c4))g (1)h 312( (c2) + (c3))g (1); h 412( (c1) + (c4))g (0)Using (3.26) the formulas (3.25) give the identical results as presented byLeijnse [1993], see (2.7){(2.8), and consequently as presented by Voss andSouza [1987], see (2.6).4. Finite volume discretizationFor a given triangulation of (polygonal) domain we construct a dualmesh of nite volumes known as Donald diagrams. In 2D case it is createdby connecting all midpoints of the edges of the element with its barycenter.By considering the integral form of (1.1){(1.2) over a nite volume Viand by applying the Green's formula we get:(4.1) RVi @t dV + R@Vi ~n ~vd = RVi QpdV ,(4.2) RVi @t( c)dV + R@Vi c~n ~v R@Vi ~n D rcd = RVi QcdVwhere ~n is a unit normal vector. Next we introduce the approximation of cand p by(4.3) ~c = NPi=1 ci(t)Ni(x); ~p = NPi=1 pi(t)Ni(x)with the piecewise linear ansatz functions on triangles and bilinear onquadrilaterals. For the volume integrals we apply a "lumped mass" ar-gument to compute(4.4) RVi @t (~c)dV @t (ci(t)) RVi dV .For the boundary integrals we apply for each line segment j @Vi bydenoting xj the midpoint of j the trapezoidal rule to compute:(4.5) Ri ~n ~vd j jj (P ci(t)Ni(xj))~n(xj) ~vj:The velocities ~vj are computed from (3.16) for the integration point in localcoordinates ( j; j):As our aim is to use nonuniform grids with possibly large elements insome parts of domain, we employ an upwind scheme for the convective termin (4.2).The discretization as brie y described here was realized in the generalsoftware toolbox UG (Unstructured Grids), see Beddies [1995] where addi-tionally the consistent velocity approximation was implemented. The result-ing nonlinear system of implicit ordinary di erential equations is then solvedby the combination of some time discretization technique and nonlinear it-erative solver, see general info on UG at http://www.ica3.uni-stuttgart.de .Although the extensive numerical experiments are only at the begin-ning, rst satisfactory results were obtained. For instance for the Elder testexample, see e.g. Oldenburg, Pruess [1995], the qualitatively stable, physi-6 cally reasonable and grid nondependent numerical results were computed.The Figure 1. presents the comparision of the consistent and noncon-sistent velocity eld that were determined from the computed concentrationand pressure distribution for t = 5 years. The gure 2. shows the concen-tration distributions for t = 4 years and 8 years.Figure 1: Consistent and nonconsistent velocity eld.Figure 2: Concentration contours.